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We show that any short-range Hamiltonian with a gap between the ground and excited states can 
be written as a sum of local operators, such that the ground state is an approximate eigenvector 
of each operator separately. We then show that the ground state of any such Hamiltonian is close 
to a generalized matrix product state. The range of the given operators needed to obtain a good 
approximation to the ground state is proportional to the square of the logarithm of the system size 
times a characteristic "factorization length". Applications to many-body quantum simulation are 
discussed. We also consider density matrices of systems at non-zero temperature. 



The application of numerical renormalization group to 
quantum systems is a natural idea with a long history. 
Despite Wilson's success with the Kondo model 1], other 
early attempts based on keeping low-lying eigenstates in 
each block were less successful?^ • The basic idea of these 
methods is to break a system into subsystems, solve each 
of the subsystems, and then join the solutions together. 

This leads to the following general question: how do 
you solve a system if you know the solution of its subsys- 
tems? Consider, for example, the following toy impurity 
problem: a single spin-1/2 impurity embedded in a band 
gap insulator. Suppose that the electrons in the insulator 
do not interact with each other, but only with the impu- 
rity spin. The problem of the electrons alone can readily 
be solved, even in an infinite system, by filling Bloch 
states up to the Fermi level, and one finds exponentially 
decaying correlations in this system. If the interaction 
with the impurity is strong, the impurity problem does 
not admit an analytic solution, but given the finite corre- 
lation length, one might approximately solve a finite pe- 
riodic region around the impurity on a computer. How, 
though, can one write down a solution for the combined 
system? How can one "sew" the two solutions together? 

In one-dimension, density matrix renormalization 
group (DMRG)j3| provides a means to do exactly this 
and has been extremely successful. The states that it 
finds are matrix product states '4','5','6'|. There exist some 
promising higher dimensional generalization of these ma- 
trix product states0, Q- Precise bounds on how well 
one can approximate a given quantum state by a matrix 
product state are still lacking, however. 

Recently, there have been several advances in under- 
standing the connection between a gap and the locality 
of correlation functions [sl 1^. ITol Ull 1 1 2j| . providing a firm 
analytical basis for the notion that a gap implies expo- 
nentially decaying correlations while a power law density 
of states implies correlations are bounded by an algebraic 
decay. In this paper, we will similarly use the existence 
of a gap to study this problem of sewing states together. 

All of the known Hamiltonians that give matrix prod- 
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uct states, such as the AKLT 13] Hamiltonian, have the 
property that the ground state is an eigenvector of each 
term in the Hamiltonian separately. We refer to this as 
a local projective Hamiltonian. The first portion of this 
paper will construct a form of arbitrary gapped Hamilto- 
nians such that the ground state is an approximate eigen- 
vector of each term separately, writing the Hamiltonian 
as a sum of terms in Eq. ijHJ). This will be a first step to 
building a matrix product form of the ground state. We 
construct such a form in this paper, but do not bound the 
number of states required in the matrix product construc- 
tion. Such a bound will be given in a future work. While 
the proofs here will be to some extent constructive, they 
will assume that certain properties of the ground state 
are known, and thus they are not so useful in themselves 
for the problem of finding ground states. 

Next, we briefiy discuss applications to numerical sim- 
ulation. Of course, one application is in analyzing exist- 
ing algorithms, but we consider suggest the possibility of 
different algorithms based on the proofs in the first part. 
In this case, we will discuss how to find the needed prop- 
erties of the ground state used in the proofs in the first 
part. 

The last portion of the paper will consider systems at 
non-zero temperature. In this case we will show that the 
density matrix of the system can be written in a matrix 
product form, which provides a higher dimensional gen- 
eralization of the one-dimensional matrix product form 
for density matrices pH- 



I. APPROXIMATE LOCAL PROJECTIVE 
FORM OF THE HAMILTONIAN 

Consider the AKLT Hamihonian|3, H = Y,H^, with 
Hi ^ ■ S^+l + il/3)iSi ■ S,+if. This Hamiltonian 
has an exact matrix product ground state. For a chain 
of A'' — 1 sites, suppose there are ground states labeled 
by an index p. Then, a chain of sites is supposed 
to have ground states labeled by an index a, with a) = 
J2i3 s ^a.i3{s)(3) (g) s), where s) denotes a complete set of 
states on the A^-th site (in this case, there are three such 
states), and Aa,p{s) is the matrix defining the matrix 
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product state. This then gives a wavefunction which sews 
the solutions of the two subsystems together. One way 
to find such wavefunctions is DMRG, while another is 
the variational matrix product method. For the AKLT 
chain, the ground state has a — 1,2 with 

For more general Hamiltonians, the range of the indices 
a, f3, s may be larger and the matrix Aa,/3{s) may be dif- 
ferent. 

This ground state not only minimizes the Hamiltonian 
H, but also minimizes each Hi individually, and thus is 
an eigenvector of each Hi. This observation is funda- 
mental to the work in this section. We take an arbitrary 
Hamiltonian and approximately rewrite it in an approxi- 
mate local projective form, defined to be a form in which 
the Hamiltonian is a sum of local terms Mi such that 
the ground state is close (as defined below) to an eigen- 
vector of each Mi separately. A Hamiltonian with such 
a form can truly be solved by solving subsystems sepa- 
rately. Break a chain of N sites up into two subchains of 
N ~m and m sites. For each chain, find the eigenvectors 
of the Hamiltonian with the correct eigenvalue (the same 
eigenvalue as the ground state of the full system has for 
the Hamiltonian of the subchain). The ground state of 
the full chain will be a linear combination of outer prod- 
ucts of the given states in each subchain. The matrix 
product method realizes this for m — 1. 

We consider an arbitrary Hamiltonian H which is as- 
sumed to have some number of degenerate ground states 
and then a gap IS.E to the rest of the spectrum. That 
is, define ^o) to be an eigenstate state of H with energy 
Ea, and let Vl/a) for a = 0...n — 1 be n distinct ground 
states while for a > n we have Ea > AE. We consider 
the case Eq = Ei = ...£"„_! = 0(141. We assume that 
H obeys the finite range conditions [slfTsll. That is, the 
Hamiltonian H can be written as a sum of terms Hi such 
that each Hi has a bounded operator norm, | \Hi\ \ < J for 
some J and such that each Hi acts only on sites within 
some interaction range R of site while there are at 
most S sites j within distance R of sites i for any i. In- 
troduce some metric on the lattice d{i,j) as the distance 
between sites i and j, while d{0,j) is defined to be the 
distance between an operator O and a site j: that is, 
the minimum, over sites i on which operator O acts, of 
d(i,j). Then, for any operator Oj which acts on a site j 
with d{i,j) > R, we have [Hi, Oj] = 0. 

We now construct the approximate local projective 
form, and in the next section construct the ground states. 
Following!^, define 

A TP l-OC 

Hf= dtm), (2) 

where 

Hiit) = H,{t) exp[-(tA£;)V(2(z)], (3) 
Hi{t) = exp(iFt)7J,exp(-ii7i), 



with q to be chosen later. That is, Hi{t) is defined fol- 
lowing the usual Heisenberg evolution of operators, while 
Hi{t) is equal to Hi{t) multiplied by a Gaussian which 
cuts off the integral in Eq.(|21) at times t of order ^/AE. 
The notation of for these operators Hi is chosen to 
indicate that we make an approximation (hence the tilde) 
to the zero frequency (hence the zero) part of the Hi . 

Now, here is the key point of the paper. We claim 
that H^, acting on a ground state, gives back another 
ground state up to some exponentially small difference. 
To see this, compute {Hf)ab = (^'a, -fff^fc), the matrix 
element of H^ between states {'^a and 5'b)_17J. A direct 
computation gives 

{H?)ab = mateM~<li^^^^rm- (4) 

Let 

Plon,= J2 
0<a<ra 

and 

Phigh^l~Plow (6) 

Thus, Plow projects onto the space of ground states while 
Phigh projects onto the remaining states. Then, from 
Eq. 101, the norm 

\PMghH^^a)\ < ||7?0||exp(-(7/2) < Jexp(-g/2), (7) 

as claimed. 

One important fact is that 

i i 

Also, Hi and iJ^ have the same matrix elements in the 

subspace of ground states: {Hi)ab = {Hi)ab if < a < n 
and {) <b < n. Finally, 

\\H^\\<m\\<j. 

The Hi are local in that the commutator of Hi with any 
operator Oj which acts only on a site j is exponentially 
small in d{i,j). This follows since Eq. (0) defines Hi as 
an integral of Hi{t) over times t; the Gaussian in Eq. Q 
cuts this integral off for sufficiently long time while for 
short time Hi{t) is local. The precise statement shown 
in the Appendix is that, for any operator Oj which acts 
only on a site j, 

m^,0,]\\ < J||O,||(g(ciZ,0 + 2exp[-(ci;Ai?)V(2g)]), 

(9) 

where I — d{Hi,j) and the function g{cil,l) is an expo- 
nentially decaying function of ?/^c for some microscopic 
length scale of order the interaction range R. The 
constant ci is a characteristic inverse velocity of propa- 
gation in the system; the existence of a finite velocity of 
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propagation, as discussed in the Appendix, is essential in 
showing that Hi{t) is local for short time. Eq. @ implies 
that H'j* is local in that it has a small commutator with 
operators which are far enough from i. 

We can further define Mi to be an approximation to 
Hi which is truly finite range: Mi will exactly commute 
with Oj if d{i,j) is greater than a certain range Iproj- To 
do this, define 

M, = J dtexph(tA£;)V2]i?f™'=(i), 

(10) 

where 

7Jf""^(t) - exp{iHi,,t)H, exp(-^i^loc^), (11) 
^loc- J2 (12) 

Thus, H\oc is the sum of terms Hj with d{i,j) less than 
Iproj — R, so that _ff*''"'^'^(t) only acts on sites within dis- 
tance Iproj of i. Thus, the procedure to define the Mi is 
very simple: one uses the definition Eq. |(2Jl, but as one 
evolves Hi(t) one drops terms which involve sites more 
than Iproj from i. 

In the Appendix, we show that 

\\H°-M,\\ < J{N{lproMci ^proj 1 ^proj ) (13) 
+ 2exp[-(ciVo,A^;)V(2g)]), 

where N{lproj) is defined to be the number of sites j with 
Iproj — R < d{i, j) < Iproj + R- Note also that | \Mi\\ < J. 

We now pick q. For a given range Iproj of the Mi, we 
want to minimize \PhighMi'i'a)\, so that the ground states 
are approximate eigenstates of the Mi. By a triangle 
inequality, 

\Ph^ghM,'fa)\ < \P|^^g|^H°^ a)\ + \\H^ - M,\\ 

< J (exp(-(7/2) + N{lproj)g{cil proj 1 ^proj ) 
+ 2exp[^{cilprojAE)y{2q)]). 

To get the best bound, we pick q = cilprojAE. Then, 

\PMghM,-^a)\ < JO{cM-lproj/lfac)), (14) 

where O denotes a quantity of order exp{—lproj /I fac), 
with Ifac being the characteristic factorization length. 
The length Ifac is equal to the minimum of (ciAE)^^ 
and ^c, and thus for small AE, Ifac = (ciA£^)~^. 
With the given q, the bound in Eq. I|13|) becomes 

— Af,;|| < J[N{lproj)g{cilproj, Iproj) 

+ 2exp{-cilprojAE/2)]. 

This difference is exponentially small in Iproj /I fac, so that 
difference between the ground state energy per site of 
H — J2i Hi and that of the Hamiltonian M = J^i ^'h 
is exponentially small in Iproj /Ifac- Defining N to be 
the number of sites i in the system, if N\\H^ — Mi\\ is 



less than of order AE', then the ground state of M has 
a non-vanishing projection onto the ground state of H. 
This requires an Iproj which is of order log(7V). 

We claim that these Mi realize the approximate local 
projective form, 

iJwA/==^M,. (15) 

i 

We start with the simplest case of only one ground state, 
n = 1. Then, Eq. (|14|l implies that the ground state "^/q) 
is close to an eigenvector of each Mi. That is, 

|Af,^'o)-(Af,;)^'o)| < JO(exp(- iproj /Ifac)), 

where (...) denotes the ground state expectation value. 
By picking Iproj large, we can make this difference as 
small as desired. 

The Mi give the desired approximate local projective 
form for the case of a unique ground state. We will 
show in the next section how to construct matrix product 
states that are approximate ground states of this Hamil- 
tonian. However, we first consider the case of multiple 
degenerate ground states. 

If there are multiple ground states with a gap to the 
rest of the spectrum, then the situation is slightly more 
complicated. Eq. j?)) implies that the Hf acting on 
ground states gives states which are close to ground 
states, but no longer necessarily implies that the ground 
states are eigenstates of the H^. Instead, it depends 
to some extent on what basis we choose for the ground 
states. As a simple example, consider the Majumdar- 
Ghosh Hamiltonian for a one-dimensional spin-1/2 chain: 

H = J:^H, with H, = JJ:,[S^ ■ S,+i + {l/2)S, ■ S,+2]- 
This Hamiltonian has two exact ground states; in one 
state sites i = 1,2 are in a singlet, sites i = 3,4 are 
in a singlet, and so on, while in the other state sites 
i = 2, 3 are in a singlet, sites i = 4, 5 are in a sin- 
glet, and so on. Denote the first state by ^even) and 
the second state by ^odd)- Now, the states 'i'even) 
and ^'odd) break translational symmetry; the expecta- 
tion value (5'cvon, ^^i^odd) is an alternating function of 
i. However, in an infinite system the expectation value 
(^odd, -ffj^'odd) vanishes. Thus, in the subspace formed 
by the two vectors 5'ovcn.odd), the Hi are diagonal and 
therefore the H^ are also diagonal in this subspace. So 
with the given Mi, there is no problem in this basis 
of ground states: the states vE'evon,odd) are approximate 
eigenstates of the Mi for large q. Of course, as is well 
known for the Majumdar-Ghosh chain, if we were to pick 
H, = {J/2)[S^ ■ S,+i + S, - S,-i + (l/2)^,_i • then 
the states ^'cvcn.odd) would be exact eigenstates of the 
Hi, but let us suppose that we do not know that this 
form of the Hamiltonian is available. 

Suppose instead we choose to form ground states which 
are eigenvectors of the translation operator by ^"5) = 
*ovc„) + *odd) and = *ovcn) - *odd)- Then, the H, 
are not diagonal in this subspace, and the states s.,a) 
are not approximate eigenvectors of the H^, no matter 
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how large q is. One way to get around this is to go to 
an enlarged unit cell of two sites, setting H — 
where flj = H2j + H-zj+i and then the states '^s,a) are 
approximate eigenvectors of Mj. However, the simplest 
solution is to use the states ^'even.odd) instead of ^'s,a)- 
Thus, the important question is: can we simultane- 
ously diagonalize all of the Mi in the subspace formed 
by the ground states '^a) for < a < n? If so, then we 
can ensure that, by the appropriate choice of basis for 
the ground states, each ground state is an approximate 
eigenvector of each of the Mi and we will have 

|(M, - < JOiexp{-lproj/lfac)), (16) 

If the states 'i'a) in this basis break translational sym- 
metry, then the expectation value of Mi in a state ^Pq) 
may depend on i even if H is translationally invariant. 
In order to simultaneously diagonalize the Mi in the sub- 
space of ground states, we need the Mi to commute in 
this subspace. 

The conditions for the Mi to commute in this subspace 
are discussed in an Appendix. We will show that, except 
for a few artificial examples, the Mi approximately com- 
mute. In particular, if H is translationally invariant, we 
will show that the commutator of the M,; in this sub- 
space vanishes exponentially in the system size, and thus 
we can pick such that Eq. H16|) holds. For H which 
are not translationally invariant, we show that most of 
the Mi commute in this subspace. More precisely, for any 
basis of the ground states, define o{i) to be the operator 
norm of the off-diagonal part of Mi in the low energy sec- 
tor in that basis (the off-diagonal part of Mi is a matrix 
which has zeros on the diagonal, but whose off-diagonal 
elements are the same as Mi). We will show in the Ap- 
pendix how choose a basis in which we can bound o(i) 
by Eq. ljB3|) . This implies that we can pick the 'i'a) such 
that 

\M,^a) - (*a,M,*a)*„)| < JO{exp{-lproj /I fac)) + o{i) , 

with the sum of o{i) bounded. 

II. MATRIX PRODUCT STATE FOR THE 
LOCAL PROJECTION HAMILTONIAN 

We pick a given ground state '^a) and approximate 
it by a matrix product state, with an appropriately 
bounded error. For simplicity, we consider the case in 
which we can pick the ^'a) such that Eq. Hlt)|) holds. This 
includes, as discussed above, all translationally invariant 
systems, as well as many translationally invariant sys- 
tems, those without local zero energy excitations as dis- 
cussed in the Appendix. We review previous work on 
matrix product or valence bond states, and then provide 
various constructions of the matrix product state in the 
given case. The value of Iproj required to obtain a good 
approximation to the ground state will be seen to grow 
as a power of the logarithm of the system size. 



A. Matrix Product and Valence Bond States 

The one-dimensional matrix product state, as dis- 
cussed above, takes a chain of — 1 sites with a set 
of ground states f3, and constructs the ground states of 
an N site chain by a) — J2p s ^a,0{s)s) x /3), where 
j3) are ground states of an Af — 1 site chain. Matrix 
product ground states in one-dimension can always be 
written as valence bond states 4] . Construct an en- 
larged Hilbert space on each site, labeling states on 
site i by two indices ai,Pi. A wavefunction is con- 
structed in the enlarged space, such that this wavefunc- 
tion '4>{ai, Pi,a2, P2, ■■■) is a product of wavefunctions 
'0(/3i, a2)V'(/32, Qfs), .... Finally, a map is written on each 
site from the original Hilbert space s) to the enlarged 
Hilbert space aipi), and the wavefunction on the orig- 
inal Hilbert space is defined to be the wavefunction of 
the mapped state on the enlarged Hilbert space. To 
make this concrete, suppose the matrix product con- 
struction gives a) — J2i3 s ^ (^)- Then, set 
■(/'(A, ttj+i) = S{Pi,ai+i). Let F map state Si) on site 
i onto J2ai0i ■ Alternately, this map F 
can be viewed as a projection from the space of states 
aif3i) onto the states s^-)|l8). Then, the matrix product 
wavefunction is given by rp(F{si, S2, ■■■))■ 

In the AKLT case, labels one of the two spin-l/2s 
and Pi labels the other one. The product wavefunction is 
a product of singlet pairs, while the map F projects the 
two spin-l/2s onto a spin-1. The matrix product state 
for the AKLT chain can be written also as a valence bond 
state. Each site has a spin-1, which may be represented 
by two spin-1/2 spins. One spin-1/2 is in a singlet with 
a spin-1/2 on the next site to the right, and one is in a 
singlet with a spin-1/2 on the next site to the left. This 
state is then projected onto the spin-1 state of the two 
spin-1/2 spins on each site. 

Matrix product states and valence bond states are 
equivalent. However, the discussion above was confined 
to pure states (wavefunctions) on finite systems. In 
matrix product and valence bond states were also con- 
structed for mixed states (density matrices) and again 
shown to be equivalent. 

Thus far we have discussed systems in one dimension. 
A higher dimensional system can always be viewed as a 
one dimensional system as follows. For a d-dimensional 
system in which each site is labeled by d coordinates, all 
of the sites with a particular value of one coordinate can 
be grouped into one supersite, leaving a one dimensional 
chain. This method is very limited in practice; for a 
system of linear size L, the size of the Hilbert space on 
a single supersite is exponential in L"^^^ and thus the 
range of the indices a, /3 is also exponentially large. This 
method amounts to studying a one-dimensional ladder 
system. 

Valence bond states are often regarded as a more ap- 
propriate way for constructing states in more than one 
dimension. To construct such a state0,0|, on an arbi- 
trary lattice in any number of dimensions, for each site i 
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one constructs one fc-dimensional auxiliary Hilbert space 
per bond, where the bond connects site i to site j. An 
wavefunction is defined in the enlarged Hilbert space, 
which is a product of wavcfunctions on each bond, where 
the wavefunction on each bond is a function of states at 
the "ends" of the bond (in one dimension, these are the 
indices Pi,ai+i). For each site, a map is defined from 
the original Hilbert space to the product of the auxiliary 
Hilbert spaces on that site. In this paper, we restrict to 
methods based on supersites for higher dimensions, while 
a future publication will provide valence bond construc- 
tions in this case 01 . 

Manipulating these higher dimensional valence bond 
states is difficult. After 4j, most higher dimensional 
work involved special examples where the Hamiltonian 
was exactly equal to a sum of projection operators 20j, 
so that the ground state is exactly a higher dimensional 
generalization of the AKLT ground state. However, in 
an important advance0, valence bond states were sug- 
gested as a good ansatz for arbitrary Hamiltonians, with 
a numerical technique being used to compute the state. 
We now provide a construction of the matrix product or 
valence bond states. 



B. Construction of Matrix Product State 

We first give the matrix product construction in one 
dimension. The approach discussed above in one dimen- 
sion is an iterative procedure: from the ground states of 
an N— 1 site chain, we construct those of an TV site chain. 
This procedure is shown schematically in Fig. 1(a). We 
first find the set of allowed states on the first two sites, 
then the set of allowed states on the first three sites, and 
so on. Our procedure will instead be a "hierarchical" 
procedure, as shown in Fig. 1(b). We will find a set of al- 
lowed states on pairs of sites; we then find sets of allowed 
states on groups of four sites, and so on. In practice in 
DMRG, the iterative approach works better, but we find 
the hierarchical approach gives better bounds here. 

Consider a one-dimensional system with given Iproj- 
First group the sites into supersites with sites I < i < 
'^Iproj grouped into one supersite, sites 2lproj + 1 < « < 
Alproj + 1 into another supersite, and so on, grouping 
2lproj sites into each supersite. Suppose each site has 
an TO-dimensional Hilbert space. The dimension of the 
space of states on the sites 1 < i < 2lproj is at most 
jj^2ipraj From now on, we refer to these supersites sim- 
ply as "sites"; with this grouping, the operators Mi act 
only on pairs of neighboring sites, with no longer range 
interaction or three site interactions. We let Mk,k+i be 
the sum of the operators Mi which act on sites k and 
k + 1. Label the supersites of the system by i = 1...N. 
The system is periodic, so that site iV -f 1 is identical to 
site i; similarly, d{N, 1) = 1, d{N, 2) = 2, and so on. 

Pick a given '^a- Let aj, label a complete basis of states 



al) on site k. From Eq. H16|l . 

|(Mfe,fe+i - («'a,Mfe,fc+i^'a))^'a)| 
^ (//pT^oj ^ (exp( ^proj / ^ f ac) ) • 

Let project onto the eigenvectors of Mk,k+i with 
eigenvalues A such that |A — {'^a, Mk,k+i'^a)\ < x, for 
some X to be chosen later. Let af. label these different 
eigenvectors af,). Then, 

X(*a,(l-Pfc)*a) < |(Affe,fe+i - (*a,Affe,fc+i*a))*a)| 
^ Jlproj ^ (^^P( ^proj / ^ f ac)) • 

Then, choose x = Jlpro]0{exp[-lproj / {'^og2{N)l fac)]), so 
that 

(*a,P,^*a) >l-6l, (17) 

where bi = exp[-(l - 1/ log2(A^))(Voi/^/ac)]- This im- 
plies that Mk^k+i has at least one eigenvalue A, such that 

|A- (^'a,Mfe,fc+l^'a)| 
< JlprojO{exp[-lproj/ {l0g2{N)lfac)])- 

Therefore, there is at least one state a|) for each site k. 

For each pair of sites k,k + I surrounded by an oval 
on the second line of Fig. 2(b), we calculate the of 
the above paragraph. Then, for each group of four sites, 
k.k + 1, k + 2, k + 3 surrounded by any oval on the third 
line of Fig. 2(b), we let P| project onto the eigenvectors 
of P^ P^_^_2Mk+l,k+2Pk+2Pk ™ t'^s space of states al) (g) 
ce'l_^_2), such that the eigenvector has eigenvalue A such 
that 

|A- («'a,Mfe+l,fc+2^'a)| 
< Jlproj[0{exp{~lproj/lfac)) + 26i] (W) . 

Let label the resulting eigenvectors a^). Define ^'^ = 
^fc^fc+2*a)/-^(*a, PkPk+2'^a)- This vcctor is normal- 
ized to unit norm. Then, {-^^-^al <2bi. Thus, 

\{Mk+l.k+2 - {'flMk+l.k+2K))'^l)\ 
^ Jlpj-QjCD (^*dxp>(^ Iproj I ^ f acj) ~t~ 2h\ Jlproj ■ 

Therefore, (*! , Plvpi) > 1 - 62, where 62 = 
(exp(-VojA/ac) + 26i)i-i/iog2W. 

Proceeding in this fashion, we find that 
(^'™-i,Pf*™-i) > 1-6™, where 6„ = 
(exp(-Voj/^/ac) + 26,„_i)i-i/iog2(Af). There are 
h = log2(A^) levels of this construction. Thus, 
after the last step, hh is bounded by a quan- 
tity of order iVexp[-(?p,.oj///ac)(l - 1/^)'"] ^ 
N exp{— Iproj / el fac)- Choose Iproj such that 

Iproj/lfac > 0(log(J/A£;)log(7V)2). Then, there is 
at least one state a^), and one may show that 

(a\ M^a'') < II + (18) 

i i 
NJO{exp[ /{\0g2{N)lfac)]). 
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a) b) 




FIG. 1: Illustration of iterative and hierarchical procedures, 
a) Iterative procedure. Filled circles represent sites, ovals sur- 
rounding sites representing grouping of sites into a supersite. 
The first stage groups sites 1 and 2. The second groups the 
combined site with site 3, and so on. b) Hierarchical proce- 
dure. The first stage groups pairs of sites, the second stages 
groups four sites into a single site, and so on. 

Thus, the energy of this state a^) will be within Ai? of 
the ground state energy. 

Thus, this procedure yields a matrix product state 
close to the ground state in energy. At the same time, 
this procedure does not yield too many states. All of the 
states a^) are within energy of order AE of the ground 
state. We can choose Iproj so that they lie within energy 
Ai?/2 of the ground state; it can then be shown that 
the number of distinct a^) is at most 2n, where n is the 
number of ground states. By choosing the bound on the 
expectation value of the energy even smaller, the ground 
state may be approximated to arbitrary accuracy. 

This procedure can be extended to systems in more 
than one dimension. One possibility is to treat a 
higher dimensional system as a one-dimensional system 
by grouping all the sites with a particular value of a co- 
ordinate into a supersite as described in the subsection 
reviewing the matrix product and valence bond construc- 
tions. Another possibility is to group group 2"^ sites to- 
gether at each stage of the hierarchy. In either case, it 
is still only necessary to take an Iproj which grows as 
a power of log(iV) to get a good approximation to the 
ground state. 

III. QUANTUM SIMULATION 

One application to quantum simulation of matrix 
product or valence bond states in higher dimensions is 
variational jS], and the results here may be useful in an- 
alyzing such algorithms. However, another application 
of these results to quantum simulation involves using the 
construction of the previous sections. The results in this 
section are illustrative, and will be worked out for more 
practical examples in a future publication^^. The goal 
here is to show, for at least some simple examples, that 
in principle one can actually use computations on finite 
systems to write down wavefunctions for much larger, or 



even infinite, systems, and to provide a variational tech- 
nique that gives a lower bound on the energy. 

To do this, we propose to calculate the matrices M^, 
determine the correct eigenvalue of each Mi, and then 
use this to determine the ground state wavefunction for 
a large system. At this point, there is a very natural ob- 
jection. The procedure requires that we find the correct 
eigenvalue of the Mi. If we can do this, then we know 
the ground state energy of the system (up to an exponen- 
tially small error). If Iproj is sufficiently large to obtain 
approximately correct eigenvalues for the M^, then why 
not just do an exact diagonalization of the system on a 
system of size Iproj and compute the ground state energy 
that way? One answer is that the ground state energy 
is not the only interesting aspect of the system. Correla- 
tion functions are much more important, and it is often 
a difficult task to determine the long-range order from 
the quantum numbers (such as momentum and spin) of 
the low-lying states found in exact diagonalization. The 
procedure we outline of sewing together solutions will 
provide a way of taking a solution on a finite size system 
and extending it to a wavefunction for a system of much 
larger, or even infinite, size. This wavefunction can then 
be used to compute long-distance behavior of correlation 
functions. Also, as we will see below, in some cases this 
yields better energy estimates than exact diagonalization 
of the same size system. 

We now discuss the first step, computing the Mi. One 
way is to directly calculate the Mi, using Eq. H1U|) and 
using Eq. Hll() to define i7|''""'^(f). One way to compute 
^trunc^^j is via a series expansion: i?f """(t) = Hi + 
it[H^°'^, Hi] + ... Another way is to exactly diagonalize a 
finite system of size Iproj , compute the matrix elements of 
Hi between the eigenvalues, and use Eq. to get matrix 
elements of iJf . In both cases it is necessary to choose q 
to minimize the difference \\H^ — Mi\\ given by Eq. Ijl^^ll 
and also the norm \PfiighHi'^a) \ of Eq. Q, giving a g of 
order cilproj^E. 

Given that the ground state is an approximate eigen- 
state of the Mi computed in this way, the question arises: 
what is the eigenvalue? That is, what is the approximate 
expectation value of (Hi) in the ground state? In some 
cases, finding the correct eigenvalue is an extremely dif- 
ficult problem. For example, consider an Ising spin-glass 
Hamiltonian: H = - Hi with Hi = JijSfSj, where 
Jij is some set of random couplings between nearby spins. 
This is a purely classical problem, since all of the Hi 
commute, and any state in which each spin has a definite 
value of is an eigenvector of every Hi. However, to find 
the ground state is clearly a difficult task! In this highly 
disordered system, each Hi has a different expectation 
value and one must find the correct eigenvalue for each 
one. For ordered systems the task is much easier. If the 
ground state does not break translational symmetry then 
each Mi has the same ground state expectation value. If 
there is a symmetry breaking ground state with an en- 
larged unit cell, there are still only a discrete number 
of different expectation values for the Mi. For example. 
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in the Majumdar-Ghosh chain, the ground states are in- 
variant under translation by two sites, and there are two 
different ground state expectation values for the Mi. 

We have performed some simple numerical experi- 
ments on systems of free particles on a lattice with dif- 
ferent gapped band structures. In this case there is no 
symmetry breaking, and if H is translationally invariant, 
then the expectation values of the Mi in the ground state 
are also translationally invariant. The lowest eigenvalue 
of Hf for given q always provides a lower bound to the 
energy as in Eq. H19|l . However, we have found that as 
q is increased, the lowest eigenvalue increased and con- 
verges rapidly to the ground state energy per site. There 
are several other eigenvalues very close in energy to the 
lowest, and then a gap to the rest, so in this case at least 
the identification of the correct eigenvalue is easy and this 
technique in fact provides a way to compute the energy 
per site. An upper bound to the energy per site as well 
as a macroscopic wavefunction can be obtained by using 
one of the matrix product constructions above. The most 
important question, of course, is how well this procedure 
can be extended to complicated interacting systems. 

Although the main point of the present paper is the 
formal construction of the Mi, we finally discuss here two 
brief attempts to apply these techniques to interacting 
systems. First, the techniques here provide a variational 
lower bound to the energy of the system. Consider a 
translationally invariant Hamiltonian, H = Hi, with 
all the Hi equal. Then, the ground state energy per site, 
Eo/N is at least equal to the smallest eigenvalue of Hi. 
Further, 

Eo/N>X^in{H,+iai[H„H]+a2m,H],H] + ...), (19) 

where Amin(O) is equal to the smallest eigenvalue of O, 
and ai, a2, ... are arbitrary constants. The operators H^ 
are given by a particular choice of the constants ak that 
can be obtained from Eq. (|2Jl. Other choice are possible, 
and varying over the constants will provide a variational 
lower bound for the energy. Thus, this provides an inter- 
esting complementary approach to other quantum simu- 
lation techniques, since almost all other techniques pro- 
vide either approximate estimates or variational upper 
bounds. 

We have studied how this bound is approached on a 
spin-1 Heisenberg chain. The ground state energy of 
this chain is known very accurately from DMRG. The 
present method H19() here is not intended to compare 
to DMRG, but rather we compare to exact diagonaliza- 
tion. In essence, this method provides another type of 
boundary condition, instead of the usual periodic or anti- 
periodic boundary condition, with the advantage that in 
this case we know rigorously how the ground state en- 
ergy compares to the energy from this procedure. We 
have some freedom to pick the constants ai,a2, .... We 
also have freedom to choose Hi to be an operator on a 
pair of neighboring sites, on three neighboring sites, or 
in general on any supercell of m sites. Since H is real 
symmetry in this case, we pick = for k odd. If 



m = 2, at the most trivial level of = for all k, we 
need to diagonalize a Hamiltonian with n — 2 sites. For 
02 7^ 0, flfe = for k > 2, we need to diagonalize a Hamil- 
tonian with n = 6 sites, and so on. If we instead we 
pick TO = 3 we need to diagonalize a Hamiltonian with 
n = 3 sites if all Ofe = 0, a Hamiltonian with n = 7 sites if 
02 7^ 0, Ofe = for k > 2, and a Hamiltonian with n = 11 
sites if 02 7^ 0, 04 7^ 0, Ofe = for k > 2. 

We have chosen to take to odd for the following reason. 
This technique provides a lower bound for the energy 
of the Hamiltonian, both for the infinite system and for 
the particular n-site Hamiltonian. However, we know for 
the spin-1 chain that on even size systems the energy 
obtained is already less than the ground state energy. 
Thus, the most efficient results will be obtained for odd 
size systems, and hence we pick an odd to. 

For TO = 3 and ak = for all k we find Eq/N > 
— 1.5, which is equal to one- half the energy of an open 
three site chain. For to = 3 and 02 ^ 0, with all higher 
Ofc = 0, we find Eq/N > —1.42569 after picking 02 — 
—0.075. This requires diagonalizing a 7 site system; the 
most naive estimate that can be obtained from a 7 site 
system, taking m — 7 gives a worse bound of Eq/N > 
-1.43909. For to = 5, 02 = -0.075 and all higher o^ = 0, 
we find Eq/N > —1.4156. Slight improvements on these 
bounds can be found by better optimization of 02. This 
last estimate requires diagonalizing a 9 site system. One 
we move to an 11 site system, we have the option of either 
considering to = 7, 02 7^ 0, or to = 5, 02 7^ 0, 04 7^ 0, with 
all higher ak vanishing. 

For purposes of computing the energy, then, this tech- 
nique offers some slight improvements over exact diag- 
onalization. Using exact diagonalization with periodic 
boundary conditions in this particular case, even size sys- 
tems offer lower bounds, while odd size systems offer up- 
per bounds, while we find a lower bound in every case. 
The estimate from a 9 site system using this method, 
for example, is better than that found using exact diag- 
onalization of an 8 site system, where the energy is esti- 
mated to be —1.417, but not as good as that found by 
diagonalizing a 10 site system, where one finds —1.4094, 
while for a 7 site system the estimate is better than that 
found from diagonalizing a six site system, where one 
finds —1.44. Further this method offers the only way to 
obtain rigorous lower bounds. This may become espe- 
cially important in studying frustrated systems. On a 
frustrated system with spiral order, for example, one has 
no foreknowledge that exact diagonalization of a peri- 
odic chain of a given length will provide a lower or upper 
estimate on the energy. The energy estimates of this un- 
frustrated chain obtained from exact diagonalization of 
odd size systems are poor compared to those obtained 
by exact diagonalization of even size systems; on a frus- 
trated chain one has no notion of which sizes will yield 
accurate estimates of the energy. 

In one dimension on unfrustrated systems, then, this 
technique does not give much improvement on the energy, 
as one can obtain a better estimate by exact diagonaliza- 



8 



tion of a system of one site more. However, in higher 
dimensions, even on an unfrustrated system, this tech- 
nique may become more useful again. Suppose we have 
an unfrustrated system of size L-hy-L in higher dimen- 
sions, and suppose it follows the pattern found here, that 
the most accurate estimate of the energy from exact di- 
agonalization is found from periodic systems with even L 
where one obtains a lower bound. The present technique 
offers the possibility of obtaining accurate estimates of 
the energy from a system of size L — 1 instead of size L, 
which means studying a system with 2L — 1 fewer sites. 

We know that the ground state wavefunction of the 
full system has a bounded projection onto states other 
than those with close to the given eigenvalue of Hi + 
a2[[Hi, H], H] + ... . Indeed, for the particular choice of 
Eq. ||2Jl, the projection onto such states is provably expo- 
nentially small in the size of the system considered. Thus, 
one can follow a procedure of breaking a chain or lattice 
into blocks, building the Mi for each block and diagonal- 
izing it in each block, restricting to the states with the 
given eigenvalues in each block, and then studying the 
behavior of the full Hamiltonian in this reduced space of 
states. As a first test of this algorithm, we take m — 3. 
If we simply take the project onto the states with low- 
est eigenvalues of the three site Hamiltonian, all = 0, 
this provides a poor approximation to the eigenstates of 
larger systems. For example, on a seven site chain, with 
the three sites taken in the middle of the chain, there 
are 3^ = 243 states such that the three site Hamilto- 
nian has an energy per bond given by E/2 — —1.5. The 
Hamiltonian of the seven site chain with open bound- 
ary conditions has a lowest energy state with energy per 
bond equal to —1.43909, but if we project onto the given 
243 states above, we only achieve an energy per bond of 
— 1.36496. However, if we take a2 = —.075 with ak — 
for k > 2, there are 243 eigenvalues of the Hamilto- 
nian with E/2 < —1.34, and then the 244-th state has 
eigenvalue —0.945706. Projecting the Hamiltonian of the 
seven site chain with open boundary onto these 243 states 
we find that an energy per bond of —1.43331. Thus, we 
have accurately selected the needed states, having pro- 
jected from 2187 states to 243 states, or from 393 states 
with Sz = to 51 states with Sz = 0. 

A further test breaking a fourteen site periodic spin- 
1 chain into seven site blocks and taking m = 3 and 
a2 = —.075 to project onto states in each block showed 
that the lowest energy wavefunction in this subspace had 
energy per site equal to —1.39545, compared to the exact 
result of —1.40394 for this size. If instead of joining sub- 
blocks we had added sites to a subblock one at a time, 
and used the present method of projecting onto the states 
in each subblock, we would arrive at an algorithm sim- 
ilar to DMRG. However, the goal is not to compare to 
DMRG in one dimension, but rather to present an algo- 
rithm that can be extended to higher dimensions, which 
DMRG cannot. 

A second technique, which may also offer the possi- 
bility of improved accuracy on the energies compared to 



exact diagonalization is based on the idea that the oper- 
ators i/f are defined by the choice of constants = 
for k odd, and = (-l)'^/2[(/j - 1)!!/A:!](g/A£;2)'= for k 
even. We then have a series expansion for in powers 
of q. We can then perturbatively expand the eigenvalues 
of Hi in powers of q and extrapolate to g = cxd from a 
finite number of terms. This is a speculative approach 
that is currently being studied. 

IV. NON-ZERO TEMPERATURE 

We now turn to systems at non-zero temperature. In 
this case the temperature enables us to construct an 
approximate matrix product form for the density ma- 
trix, regardless of whether or not there is a gap. The 
unnormalized density matrix for the system is equal to 
p — exp{~f3H). 

We will construct an approximate matrix product 
form, p{l3,lproj), so that 

p{f3) « p(/?, ) (20) 

= ^ Pi(ai)p2(a2)...Fi({a,})i^2({a,})... 

{"^ 

where for each site i we assign an index ai defined be- 
low, and sum over all values of that index. The operator 
Pi{ai) acts only on site i, and the functions Fi obey a 
finite range constraint: each Fi depends only on the aj 
for d(i,j) < Iproj for some Iproj defined below. The error 
between p(/9) and p{l3, Iproj) will be exponentially small 
in Iproj , while the range of the indices Oii will depend on 
Iproj- Specifically, we bound the error by showing that 
for any operator O, Z-i|Tr[Op(/3)] - Tr[Op(/3, Voj)]| < 
c||0||, for some constant c, where we define Z = Tr(p(/3)). 

Before defining p{P, Iproj), we recall the Trotter- 
Suzuki[23| decomposition of the path integral. In 
this case, we write « p„, where p„ = 

[Hi exp(— /Ji/i/n)]", where the product ranges over all 
i in some given sequence; since the different Hi do not 
commute, the result depends on the particular sequence 
chosen. We claim that each of the p„ can be written 
exactly in a matrix product form as in the right-hand 
side of Eq. (IJOJ- A given operator exp(— /3i?i/n) acts on 
sites within a distance R of site i, and can be written as 
exp(-/377,/n) = E{a..,} nj,d(ij)<fl 
where the operator Oj acts only on site j and the range of 
values of index aij is exponentially large in S, the num- 
ber of sites within distance R of site i. Here, Fi{{ai,j}) 
is some function of the S different indices ai,j with the 
given i. The operator p„ is a finite product of these op- 
erators exp(— /3i/i/n). Each term in this product can be 
written in matrix product form. For each site i, the oper- 
ator eKjp{—(3Hi/n) appears n times in the given product, 
and on the m-th time it appears, we use a set of in- 
dices to provide the matrix product form as above. 
For each site, j, we have n different indices a™- for each 
site i within distance R of site j, and thus at most nS 
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FIG. 2: Example of a set of active bonds shown as solid lines; 
dashed lines represent bonds which are not active while circles 
represent sites. The Hamiltonian is a sum of terms which act 
on pairs of neighboring sites, so that the bonds connect only 
neighboring sites. Using a Manhattan metric for the lattice, 
the set of active bonds shown here is a term in p(/3, Iproj ~ 2) . 
There are four distinct clusters, three with diameter two and 
one with diameter one. 



different indices for each site j. Grouping all such in- 
dices for a given site j into one index aj, and defin- 
ing Fi{{a}) = rim^i Piii^^i^j})! arrive at the matrix 
product form for p„. 

Thus, the Trotter decomposition gives an approximate 
matrix form, but the error in this approximation com- 
pared to the exact result is not very good. In contrast, 
the stochastic series expansion|23 provides a much better 
way of approximating the desired pi/3) with much smaller 
error, but it is difficult to write the stochastic series ex- 
pansion result in a matrix product form. Below, we will 
propose a matrix product form with a bound on error 
comparable to that in the stochastic series expansion. 



A. Percolation Transition at High Temperature 

The exponential exp(— /3i7) can be expanded as a 
power series 1 — (3H + ... We will show that at suffi- 
ciently high temperatures, f3~^ ~ J, there is a "percola- 
tion transition" in this exponential, as we now describe. 
Any given term in the power series expansion is a prod- 
uct of Hi for different sites i. For each term, we define a 
set of "active bonds" and "clusters" as follows: for each 
Hi which appears in the given term, we connect by ac- 
tive bonds all sites acted on by Hi, so that the length 
of the active bonds is at most 2R. We define a cluster 
to be a set of sites, all connected to each other by ac- 
tive bonds, and not connected to any other sites outside 
the cluster by active bonds. Then, for P < (3o, where (3o 
is specified below, define piP, Iproj) to include only the 
terms in the power series such that no two sites with 
d{i,j) > Iproj , are in the same cluster. Thus, each term in 
p{l3, Iproj) is a product of operators, each operator acting 
on the sites within a given cluster, such that each cluster 
has a diameter at most Iproj- See Fig. 2. For (3 < /3o, at 
temperatures above the percolation transition, we will be 
able to bound the difference between p{/3) and Iproj)- 

Let C be some set of sites i. Define He = J2i£C 



Define B{C) to be equal to the set of all sites j such that 
j ^ i for any i £ C and such that Hj acts on a site i for 
some site i £ C. Then, H^i^Q) — J2i£B(c) Then, if a 
given term in the power series expansion for exp(— 
includes Hi for every i G C, but does not include Hj for 
all j G B{C), then this term in the power series expansion 
has a cluster which includes exactly the sites acted on by 
Hi for all i € C. 

For any operators Oi, O2, 

Z-iTr(OOi(iri)02(jT2)...exp[-/3i/]) < ||0||||0i||||02||..., 

where 0{iT) = exp{—HT)Oexp{HT) and < ti < r2 < 
... < p. Thus, if C has nc sites and B{C) has nB{c) 
sites, we find that for any operator O, 

Z-'Tv{Oexp[-PiH-Hc-HBic))]) (21) 
< \\0\\exp[{nc + nB{c))jP], 

as can be seen by using a power series expansion of 

Z-iTr(Oexp[-/3(ff ~ He ~ ^5(0)]) 

rP 

= Z-iTr(Or exp[ / AT{Hc{ir) + HBi^c){^r))\ exp[-/3i7]), 
Jo 

where T gives the r-ordered exponential. 

We now define p(J3, C) to be equal to the sum of all 
terms in the power series expansion of exp(— /3_ff) that 
include Hi for all i e C, but do not include Hj for any 
3 e B{C). For any O, 

Tr(Op(/3,C)) < [exp(J/3) - 1]"^ x 
exp[/3| |iJc + -ffB(c) I |]Tr(0 exp[-/3i/] , 

as can be seen by a power series expansion. Using | \Hc + 
tlB(C)\ \ < {nc +nB{C})J < SJnc gives 



(22) 



Z-'Tt[Op{P,C)] 
< ||0||[exp(J/3) - 1]"^ exp[ncSjp]. 

Then, p{P) — p(/3, Iproj) is a sum over terms which in- 
clude a sequence of active bonds connecting any two sites 
i,j separated by a distance at least I 



proj - 



This differ- 
ence p{P) - p{P, Iproj) = - Z]m=l(-l)'"Pm> whcrC p„i is 

equal to the sum of all terms in the power series expan- 
sion of p which include at least m clusters with diam- 
eter greater than Iproj- Using Eq. (j^ : Z~^Tv{Opi) < 
X^ciustors y""^ 11*^1 1 J where the sum ranges over all different 
clusters on the lattice which connect two sites separated 
by distance greater than Iproj and nc is the number of 
sites in the cluster and where y = [exp( J/3) — 1] exp[S'J/3]. 
Similarly, Z-^Tr{Op„,) < (l/m!)(Eciustersy"^)"l|0||. 

We now show the "percolation transition" . The num- 
ber of clusters of nc sites is bounded on a regular lattice 
by Nx"'^ for some constant x, as is known for lattice 
animals 24] , where N is the number of sites on the lat- 
tice. Thus for sufficiently small 

00 00 

(1M!)( E y"'')" ^ exp[7V E {xy)-] - 1. 



m—1 



clusters 
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For fixed J,S, we can make y as close to zero as desired 
by taking suSiciently small (3. Thus, there exists some 
small enough (3o such that for (3 < f3o, 

oc 

^(l/m!)( Y y"^r<exp[7Vexp(-Zp,„,/eperc)]-l, 

m— 1 clusters 

where ^perc is of order the interaction range R. Then, 

Z-'\Tv[Opm - Tv[Op{PJproj)]\ (23) 

< {exp[N exp{-lpro] / £.perc)] - 1}||0||. 

B. Matrix Product Form at High Temperature 

We now show that p{P, Iproj) realizes the matrix prod- 
uct form Eq. H20|) for P < Pq. Each term in p{P, Iproj) is 
a product of terms acting on clusters of sites. We now 
specify the indices at. First, if site i is in a cluster, the 
index indicates all the sites in that cluster. There is 
some redundancy in this description: if two sites i and j 
are in the same cluster, then the indices and aj spec- 
ify the same set of sites. Each will also specify some 
additional information. The sum of all terms which in- 
clude a given cluster is some operator acting on the sites 
in the cluster; this operator can be decomposed into a 
sum of products of operators that act on individual sites 
in the cluster. The indices will also keep track of the 
terms in this sum. The number of different values that 
each index ai can assume is exponentially large in 1^^.^^ in 
this construction. However, this is much smaller than the 
number of different values that would be needed to de- 
scribe a general operator p: in general, one might need a 
range of values which is exponentially large in the system 
size. Thus, p(/3, Iproj) can be exactly written in a matrix 
product form, as in the right-hand side of Eq. H20I) . 

Further, we have shown that p(P, Iproj) is a good ap- 
proximation to p{P). This can be expressed in terms of 
the trace norm. For any operator O, the trace norm \0\ is 
equal to Tr(vO^O), where Tr(...) denotes the trace over 
all states and the unique positive square-root is taken. 
For a Hermitian operator, such as p{P) or Iproj), \0\ 
is equal to the sum of the absolute values of the eigenval- 
ues. For a positive definite Hermitian operator, such as 
p{P), the trace norm is equal to the trace. Use Eq. if^ . 
valid for arbitrary O, and choose O so that in a basis of 
eigenvectors of p(/3) — p(/3, Iproj), O is a diagonal matrix 
with each diagonal entry equal to ±1, the sign being cho- 
sen the same as the sign of the corresponding eigenvalue. 
Then, 

\p{P)-p{P,lpro])\ < {exp[N exp{-lproj /^perc)]-l}\p{P)\■ 
{24) 

By taking an Iproj that is of order \og{N)^perc, we can 
obtain a small error |p(/3) — p{P, lproj)\- 

We now obtain a matrix product form at arbitrary 
temperature /3 that guarantees positive definiteness of 
p{P, Iproj)- For P > Po, set 

p{P, Iproj) = piP/nJproj)"', (25) 



where n is the smallest even integer such that P/n < Pq. 
By picking n even we guarantee that p{P, Iproj) is positive 
definite. Eq. (|25|l realizes a matrix product form l|20|l : 
for each site i there are n different indices ai, one from 
each of the p{P/n, Iproj)- We group these into one single 
index, and call the new index a;, giving Eq. (|20|l . The 
number of different values that the new ai can assume 
is exponentially large in {P/Po)iproj- We can bound the 
error at arbitrary temperature as follows. 

For each p{P/n, Iproj) we define clusters. We can show, 
as above, that 

Z-^TT{Op{kp/n)p{p/n, C)p{{n -k- l)p/n)) < \\0\\y'"' 

for any < k < n — 1, where y = [exp(J/3/n) — 
l]exp[Sjp/n]. 

Define pm to be the sum over all terms in which there 
are at least m clusters with diameter greater than Iproj- 
Then, 

z-iTr(Op„o < E "y"'^rilo||- 

clusters 

Following the same arguments as before 

Z-'\Tt[Op{P)] - Tr[Op(/3, Vo,)]| (26) 
< {exp[(/3//3o)iVexp( /^perc)]~l}\\0\\. 

and so we need an Iproj that grows logarithmically in 
the system size to obtain a good approximation. At the 
value of Iproj that gives a good approximation, Iproj ~ 
X) 'Cperc log(iV), the number of different values that each 
index a can assume is of order 

exp[Oi^^,rMN)''iP/Po)]. (27) 

The final matrix product form turns a quantum sta- 
tistical mechanics problem into a classical statistical me- 
chanics problem, as tracing over the pi gives a probability 
distribution for the ai. However, the resulting classical 
problem may suffer from a sign problem. 

V. DISCUSSION 

We have shown how to construct an approximate lo- 
cal projective form of short-range, gapped Hamiltonians, 
and used this to show that the ground states are close 
to matrix product states. The main goal of the present 
paper is at a formal level; the techniques developed in 
this paper provide a way to, in principle at least, extend 
calculations on small systems to wavefunctions on much 
larger system sizes. It is worth comparing to '25^, which 
showed how to write a very general class of Hamiltonians 
as a sum of projection operators. The projection opera- 
tors there were two-by-two matrices which makes it much 
easier to find the ground state. Here, the projection op- 
erators are large matrices, and the task of constructing 
the ground state from these operators is tricky. However, 
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the important advance here is that the projection opera- 
tors are local. This strongly constrains the ground state 
and leads to the matrix product for the ground state. 

In addition to the formal interest, this work may find 
practical use in numerical simulation, as, at least in sim- 
ple cases, it is possible to directly calculate the local pro- 
jective form. Further, we have provided (|19|l . a varia- 
tional lower bound on energy. 

We finally return to the impurity problem raised in 
the introduction. The following procedure will generate 
a good wavefunction for the whole system. First, nu- 
merically compute the Mi for some region around the 
impurity. Next, compute the appropriate Mi for the sys- 
tem outside that region; this can be done analytically 
since that part of the system is non-interacting. Then, 
determine the appropriate eigenvalues for the and get 
a basis of states aimpurity) 'S) ainsuiator)- Finally, use the 
Mi that connect the two regions to determine a wave- 
function in this basis of states. This procedure can be 
followed even if there are many impurities embedded in 
the system and the resulting wavefunction can be used 
as a starting point for further improvement. 

Acknowledgments — This work was supported by DOE 
contract W-7405-ENG-36. 



APPENDIX A: LOCALITY 

The locality of the relies on the finite group ve- 
locity result 'g, [ij, 1231 . Given the finit e-range conditions 
on the Hamiltonian above, one can bound the commuta- 
tor ||[A(t),B(0)]||, where A{t) = exp(im)Aexp(-im), 
and show that this commutator is exponentially small for 
times t less than cil where I is the distance between A 
and B and ci is some characteristic inverse velocity. The 
bound is that \\[A[t),Bm\ < \\A\\\\B\\j:, 9{t,d{A,j)), 
where the sum ranges over sites j which appear in op- 
erator B and where the function g has the property 
that for |t| < cil, g{cil,l) is exponentially decaying in 
I for large I with decay length ^p. Also, for t < t', 
git,l)<it/t')git',l). 

Consider [H^,Oj] for some Oj which 
acts only on site j. This equals 

{AE/V2^) dm{t),0,] exp[-(tAi?)V(2g)]. 
Applying a triangle inequality to the integral we have 



2J\\Oj\\. Thus 



\\[H°,OM < {AE/^q) / At\\mt),OM\ X 

J — OO 

exp[-(iA£;)V(2g)]. 

Let I = d{Hi,j) be the distance between Hi and site 
j. We split the integral over times t into a sum of one 
integral over \t\ < cil and one integral over \t\ > cil. For 
\t\ < cil, we use the finite group velocity bound to bound 
\\[H,{t),Oj]\\ while for \t\ > al we use \\[H,{t),Oj]\\ < 




\\[HlO,]\\<m 
gicil,l)eM-(tAE)y{2q)] + 

exp[-(iAi?)V(2g)]) < 



/^7rq j\t\>cii 

J||O,||(5(ci/,0 + 2exph(ci;A£;)V(2<z)]), 

giving the claimed bound on the commutator 0. 
We now bound the difference — Mi\\. We have 

I - M,| I < {AE/y/2^) X (A2) 

J dtexp[~{tAEf/{2q)]m{t) - 

The difference between Hi{t) and is due to the 

different Hamiltonians used to define the time evolution. 

We now bound the difference between Hi{t) and 
^trunc^^^ This result will also be needed in the non- 
zero temperature calculation. We can replace Eq. H12() 
by iJioc = Y.j,d{i,j)<ip^,,-R.Hj + Lj,d(ij)>v„j+fl^i; by 
adding the terms Hj with d(i,j) > Iproj -I- i? we do not 
change iJf""'=(t). Then, 



\mt)-Hr^\t)\\ 
< f dt'\\[H,,Hm\- 

3.lproj-R<d{i,3)<lj,ro3+R 

Using the finite group velocity bound we have 

2R), 

j,lproj—R<d(i,j<lproj+R ^ 

where Iproj — 2R is the minimum distance between the 
operators Hj and Hi and S is defined to be the number 
of sites acted on by the operator Hj. However, the inte- 
gral JS dt'g{t', Iproj — 2R) is bounded by g{t, Iproj)^- 
Thus, 

\mt)~Hr^^^m< E jgit.iprom 

j,lpraj—R<d{i,j)<lproj+R 

— ^ {Iproj^ J g(i 7 lproj\ 

where N {Iproj) is defined to be the number the number 

of sites j with Iproj ^ R < d{i,j) < Iproj + R- 

Eq. (|A3I) can be expressed in a more general form. 
For any operator Oi, we can define Of'""'^(t) by 
Qtrunc ^ exp(^i^loc^)0^exp(-^i^loc^), with iJioc = 
Sxd(^J)<^*™„.-i^,^i• Then, Or'"'%t) only involves sites 
within distance hrunc of i (assuming that Oi originally 
involved only sites within that distance). Then, 

\mt)-or^'^m (A4) 

E \\0\\g{t,ltrunc) 

-R.<d(i,j)<ltr^^r_+R 

= N{ltrunc)\\0\\g{t,ltrunc), 



< 



3,h 
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where N(ltrunc) is defined to be the number the number 

of sites j with kmnc - R< d{i,j) < krunc + R- 

In Eq. IIA2|l . for \t\ < cilproj we can use the bound on 
g to bound the difference \ \H^{t) - i7f |, while for 
\t\> cilproj we use \ \H^{t) - < 2 J. Thus, 

11^° - M,\ \ < JiN{lproj)9{cilprojJproj) + (A5) 

2exp[-{cilprojAE)y{2q)]), 

giving Eq. ((1^ . 

APPENDIX B: COMMUTATOR IN LOW 
ENERGY SECTOR 

Here we consider sufficient conditions for the Mi to 
commute in the low energy sector. We consider the case 
that n is uniformly bounded above, and AE is uniformly 
bounded below, independent of system size. The impor- 
tant result is that for most systems of interest, including 
all translationally invariant systems, the commutator of 
the Mi is exponentially small in the system size; for a few 
examples certain of the Mi do not commute, but even in 
that case the lack of commutation only poses problems 
in a finite region of the system (in this case, as we will 
see, there are local zero energy degrees of freedom). 

The operators Mi are local, in that each acts only 
within a finite range Iproj of site i. Consider then the 
operators M°, for some q. This q is not necessarily the q 
of Eq. Using the locality results, we have 

||[M°,M,]|| < J^(NMg{cil,l) + 2exp[-{cilAE)yi2q)]), 

(Bl) 

where I — d{i,j) — 2lpj.oj is the smallest distance between 
sites acted on by Mi and Alj and Nm is the total number 
of sites acted on by Mj. 

Now, PlowM^Plon, = PlowM.Pio^. 

Therefore, [Plow Plow. Plow MjPi„J\ 

[Plow Mi Plow, Plow Mj Plow]- The second commuta- 
tor is the low energy commutator that we wish to 
evaluate. We can bound the first commutator by a 
triangle inequality, 

[[[PiowM° Plow, Plow M,Piow]\\ < ||[Af°,M,]|| 

+ {[[PMghM?Plow[[ + [[PlowMl'PMgh[[)[[M,\\. 

Following Eq. O, [[P^ghM^PioaM < Jexp(-g/2) and 
also [[PiowM^PhighM < Jcxp(— (7/2). Combining this 
with Eq. (|BT|| . 

[[[PiowM,Piow,PiowMjPiow]\\ < j'^{NMg{cil,l) 
+2 exv[-{cilAEf/{2q)] + 2exp[~~q/2]). 

Pick q = cilAE, we find that the commutator of 

1 1 [PiowM.Piow , PiowM.PlOw] 1 1 (B2) 
< J^{Nm9{ciI,1) + Aeyip[-cilAE/2]). 



Thus, the commutator of Mi with Mj in the low energy 
sector is exponentially small in d{i,j). 

For a translationally invariant, c?-dimensional, system 
on a lattice, this suffices to show that the commutator 
of Mi with Mj in the low energy sector is exponentially 
small in the system size. We can introduce coordinates 
for each site: {x,y). The c?-dimensional vector x labels 
the particular unit cell to which the site belongs, while 
the coordinate y indicates the particular site within that 
cell. The Hamiltonian is invariant under translation so 
we can pick a basis in which the n ground states are 
eigenvectors of the translation operators. In this basis, 
the matrix element elements of the M are related by 
{Ms^y)ab = {Ms+s,^y)abexp[ix' ■ {ka - h)], where ka is 
the momentum of ground state a, for < a < n — 1. 

In the simplest case (a case which applies to every sys- 
tem of which we are aware) , the ka are vectors of rational 
multiples p/q of 27r, such that p, q are integers with some 
denominator q which is independent of system size for 
large enough systems. That is, 4 = 27r(p^/g^, ...,pi/q^). 
Then, it is possible to identify a supercell, such that all 
of the n ground states are unchanged by translation by a 
supercell. The size of this supercell in a given direction 
m is equal to the least common multiple of the n dif- 
ferent g™. For example, in the Majumdar- Ghosh chain, 
the two lowest states have momentum and tt, so that 
in one case p — 0, q — 1, and in the other p — l,q — 2. 
Then, the ground states are unchanged under translation 
by two sites. For a system with momenta 0, 7r/3, 27r/3 for 
the ground states, the ground states would be unchanged 
under translation by three sites. Then, consider any com- 
mutator [Msi,yi,Ms^^y^]. This equal [Ms^,yj^, Mg^+gi^y^] 
where x' translates X2 by some number of supercells. 
By choosing x' correctly, we can make the distance be- 
tween (afi,?/2) and {x2 + x',y2) of order the linear size 
of the system. Then, we can use the result that the 
commutator vanishes exponentially in the spacing be- 
tween the operators to show that the original commuta- 
tor [Mg, i„ , My „ J is exponentially small in the system 
size, as desired |27|. 

A slightly more complicated, but very artificial, case, 
is that in which the p, q depend on system size in such a 
way that the size of the supercell is equal to the linear 
size of the system. It is not clear that such a thing can 
actually happen in a system with a finite number n of low 
energy states. Even in this case it is possible to show the 
commutator of any two Mi, Mj in the low energy sector 
is exponentially small in the system size, by using the 
smallness of the commutator [Mi,Mj] for large d{i,j) 
and expanding the commutator in intermediate states. 
However, the proof is sufficiently artificial that we do not 
give it here. 

What if a system does not have translational invari- 
ance? In this case, the Mj need not commute in the 
low energy sector. Consider the following example sys- 
tem, a one-dimensional system of spin-1/2 on each site i. 
For i ^0,1, 2, we have Hi = Sf , while Hq = + S^, 
Hi = -SS - Sf, and H2 = -S^ + S^. This is a com- 
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plicated way of writing the Hamiltonian H = X^i^^o ^ 
since all the terms acting on site cancel between iJo, 
and H2- This Hamiltonian has a doubly degener- 
ate ground state, with all spins pointing down, except 
for spin which can point in either direction. For large 
enough q, Iproj, the commutator of the Mi in the low en- 
ergy sector is close to that of the Hi in the low energy 
sector. However, clearly Hq and Hi do not commute in 
the low energy sector, since [Sq,Sq] ^ 0. Similarly, Hi 
and i?2 do not commute, and Hq and H2 do not com- 
mute. Still Hq commutes in the low energy sector with 
Hi for i > 2. 

In this case, the lack of commutation is localized near 
site 0. The two ground states differ only locally, on site 
0, and only the Hi for i — 0, 1, 2 fail to commute, while 
the others are diagonal in the low energy subspace. 

We now show how to choose a basis of the ground 
states to bound o(i). Consider a system with n — 2. 
Then, find the j which maximizes SXj, where SXj is 
the difference between the two eigenvalues of Mj in 



the low energy sector, and work in a basis which di- 
agonalizes Mj in the low energy sector. Then, o{i) < 
\\[PiowM^Pio^,Pio^MjPioyj]\\/iSXj). FromEq. the 
commutator is an exponentially decaying function of 
d{i,j), and summing over all i we find that 

o{i) < 0{J^/cilAE5\j), (B3) 

i 

independent of system size. 

For a system with n > 2, we can proceed similarly; 
we first find the Mj with the maximum difference be- 
tween its largest and smallest eigenvalues, and and use 
this to bound the off-diagonal matrix elements of other 
Mi between the corresponding eigenvectors. We then 
find an M^ which has two different eigenvalues, with dif- 
ferent eigenvectors from Mj , and show that off-diagonal 
matrix elements between those two states are exponen- 
tially small in c?(i, k). We proceed like this until we have 
bounded all off-diagonal matrix elements. 
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